Landau distribution

Landau distribution#

The Landau distribution is a highly right-skewed, heavy-tailed continuous distribution that appears famously in high-energy physics: it models the stochastic energy loss of a fast charged particle traversing a thin absorber (“Landau straggling”).

Unlike many familiar distributions, the Landau distribution has no finite mean or variance, so estimation and testing should lean on robust summaries like the median, quantiles, and likelihood-based methods.

Learning goals#

  • Know the definition (PDF as an integral) and the support/parameters.

  • Build intuition from the energy-loss origin and connect it to stable laws.

  • Understand which moments do not exist and which summaries still do (median/quantiles).

  • Implement NumPy-only sampling via the Chambers–Mallows–Stuck method.

  • Visualize the PDF/CDF and typical Monte Carlo behavior.

  • Use scipy.stats.landau for pdf, cdf, rvs, and fit.

import platform

import numpy as np

import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy import optimize
from scipy.stats import landau, moyal, norm

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(7)

print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)

Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0

1) Title & classification#

  • Name: landau

  • Type: continuous

  • Support: \(x \in (-\infty, \infty)\)

  • Parameter space (SciPy): loc \(\in \mathbb{R}\), scale \(>0\) (no shape parameters)

SciPy’s parameterization uses the usual location/scale transform:

\[f(x;\,\mathrm{loc},\mathrm{scale}) = \frac{1}{\mathrm{scale}}\,f\!\left(\frac{x-\mathrm{loc}}{\mathrm{scale}}\right).\]

2) Intuition & motivation#

What this distribution models#

Landau derived this distribution (1944) as an approximation to the energy loss by ionization of a fast charged particle in a thin layer of material.

  • The particle undergoes many small energy transfers (ionizations).

  • Occasionally, it produces a rare large transfer (a “delta ray”).

  • Those rare events create a long right tail: extreme losses happen infrequently but are much larger than typical losses.

Typical real-world use cases#

  • High-energy physics: modeling \(\mathrm{d}E/\mathrm{d}x\) in tracking detectors.

  • Radiation/particle instrumentation: charge deposition in silicon sensors.

  • Any setting where outcomes are mostly moderate but have rare, very large positive excursions.

Relations to other distributions#

  • Stable laws: Landau is a special case of a \(\alpha\)-stable distribution with stability index \(\alpha=1\) and maximal right skew.

  • Cauchy: also \(\alpha=1\) stable, but symmetric (no skew). Landau is its strongly skewed cousin.

  • Moyal distribution: a convenient analytic approximation to Landau often used in detector physics.

3) Formal definition#

PDF#

A standard (“unit”) Landau random variable has density

\[f(x) = \frac{1}{\pi}\int_0^\infty \exp\bigl(-t\log t - x t\bigr)\,\sin(\pi t)\,dt, \qquad x\in\mathbb{R}.\]

This integral definition is typical: the Landau PDF has no simple closed form in elementary functions.

CDF#

The CDF is defined in the usual way:

\[F(x) = \int_{-\infty}^{x} f(u)\,du.\]

In practice, \(F\) and \(f\) are evaluated numerically; SciPy delegates these computations to the Boost special functions implementation.

4) Moments & properties#

Mean/variance/skewness/kurtosis#

  • Mean: does not exist (diverges).

  • Variance: does not exist (diverges).

  • Skewness/kurtosis: undefined because they require finite moments.

Robust summaries that do exist and are useful:

  • Median and other quantiles.

  • Mode (maximum of the PDF).

MGF / characteristic function#

  • The MGF \(M(t)=\mathbb{E}[e^{tX}]\) does not exist for any \(t>0\) due to the heavy right tail (so it is not defined in a neighborhood of 0).

  • The characteristic function \(\varphi(t)=\mathbb{E}[e^{itX}]\) exists for all \(t\).

A standard Landau characteristic function can be written (with \(\varphi(0)=1\)):

\[\varphi(t)=\exp\left(-|t|\left[1 + i\,\frac{2}{\pi}\,\operatorname{sign}(t)\,\log|t|\right]\right).\]

Entropy#

There is no commonly used simple closed form; SciPy can evaluate it numerically.

# Basic numerical summaries for the *standard* Landau (loc=0, scale=1)

mean, var, skew, kurt = landau.stats(moments="mvsk")
median0 = float(landau.median())
q25_0, q75_0 = landau.ppf([0.25, 0.75])
iqr0 = float(q75_0 - q25_0)
f0_at_median = float(landau.pdf(median0))
entropy0 = float(landau.entropy())

print("mean, var, skew, kurt:", mean, var, skew, kurt)
print("median:", median0)
print("IQR:", iqr0)
print("pdf(median):", f0_at_median)
print("entropy:", entropy0)

# Mode via 1D optimization of the logpdf
res = optimize.minimize_scalar(lambda x: -landau.logpdf(x), bracket=(-5, -0.5, 5), method="Brent")
mode0 = float(res.x)
print("mode (approx):", mode0)

mean, var, skew, kurt: nan nan nan nan
median: 0.5756301439450783
IQR: 2.9685804468711794
pdf(median): 0.20448174823121082
entropy: 2.3726364400044817
mode (approx): -0.42931451047613484

def landau_cf_standard(t: np.ndarray) -> np.ndarray:
    """Characteristic function of the standard Landau distribution."""

    t = np.asarray(t, dtype=float)
    out = np.empty_like(t, dtype=np.complex128)
    mask0 = t == 0
    out[mask0] = 1.0 + 0.0j
    tt = t[~mask0]
    out[~mask0] = np.exp(
        -np.abs(tt)
        * (1.0 + 1j * (2.0 / np.pi) * np.sign(tt) * np.log(np.abs(tt)))
    )
    return out


# Visualize Re/Im of φ(t)
t = np.linspace(-12, 12, 3001)
phi = landau_cf_standard(t)

fig = make_subplots(rows=1, cols=2, subplot_titles=("Re φ(t)", "Im φ(t)"))
fig.add_trace(go.Scatter(x=t, y=np.real(phi), mode="lines"), row=1, col=1)
fig.add_trace(go.Scatter(x=t, y=np.imag(phi), mode="lines"), row=1, col=2)
fig.update_xaxes(title_text="t", row=1, col=1)
fig.update_xaxes(title_text="t", row=1, col=2)
fig.update_layout(width=950, height=350, showlegend=False)
fig.show()

print("|φ(t)| = exp(-|t|)  (the magnitude ignores the log phase term)")

|φ(t)| = exp(-|t|)  (the magnitude ignores the log phase term)

5) Parameter interpretation#

SciPy uses the standard location/scale transform:

\[X \sim \mathrm{Landau}(\mathrm{loc},\mathrm{scale}) \quad\Longleftrightarrow\quad X = \mathrm{loc} + \mathrm{scale}\,Z, \; Z\sim\mathrm{Landau}(0,1).\]

So:

  • loc shifts the distribution horizontally.

  • scale stretches it: quantiles, median and mode all scale linearly.

A common alternative parameterization#

Because the Landau law is 1-stable, many references use a parameter pair \((\mu, c)\) in which the stable addition rules are simple. In that parameterization, the scale also induces a location correction when converted to SciPy’s (loc, scale):

\[\mathrm{loc} = \mu + \frac{2c}{\pi}\log c, \qquad \mathrm{scale}=c,\]

or equivalently

\[\mu = \mathrm{loc} - \frac{2\,\mathrm{scale}}{\pi}\log(\mathrm{scale}).\]

We’ll exploit this in the generative section when adding independent Landau variables.

# How loc and scale change the shape


def plot_landau_parameter_effects() -> None:
    x = np.linspace(-5, 20, 4000)

    fig = make_subplots(
        rows=1,
        cols=2,
        subplot_titles=("Varying loc (scale=1)", "Varying scale (loc=0)"),
    )

    for loc in [-2.0, 0.0, 2.0]:
        fig.add_trace(
            go.Scatter(x=x, y=landau.pdf(x, loc=loc, scale=1.0), mode="lines", name=f"loc={loc}"),
            row=1,
            col=1,
        )

    for scale in [0.5, 1.0, 2.0]:
        fig.add_trace(
            go.Scatter(x=x, y=landau.pdf(x, loc=0.0, scale=scale), mode="lines", name=f"scale={scale}"),
            row=1,
            col=2,
        )

    fig.update_xaxes(title_text="x", row=1, col=1)
    fig.update_yaxes(title_text="density", row=1, col=1)
    fig.update_xaxes(title_text="x", row=1, col=2)
    fig.update_yaxes(title_text="density", row=1, col=2)
    fig.update_layout(width=1050, height=420)
    fig.show()


plot_landau_parameter_effects()

6) Derivations#

Tail behavior (key fact)#

A crucial asymptotic for the standard Landau density is the power-law right tail:

\[f(x) \sim \frac{2}{\pi x^2} \quad (x\to\infty), \qquad \mathbb{P}(X>x) \sim \frac{2}{\pi x}.\]

This is the source of divergent moments.

Expectation (why the mean does not exist)#

A (finite) mean requires

\[\int_{0}^{\infty} x\,f(x)\,dx < \infty.\]

But with the tail approximation \(x f(x) \sim \tfrac{2}{\pi x}\), we get

\[\int^\infty \frac{1}{x}\,dx = \infty,\]

so the mean diverges (logarithmically).

Variance (why it does not exist)#

Similarly, \(\mathbb{E}[X^2]\) would require

\[\int_0^{\infty} x^2 f(x)\,dx < \infty,\]

but \(x^2 f(x) \to 2/\pi\), so the integral diverges like \(\int^\infty 1\,dx\).

Likelihood (loc/scale)#

For i.i.d. data \(x_1,\dots,x_n\) from \(\mathrm{Landau}(\mathrm{loc},\mathrm{scale})\), the log-likelihood is

\[\ell(\mathrm{loc},\mathrm{scale}) = \sum_{i=1}^n \log f\!\left(\frac{x_i-\mathrm{loc}}{\mathrm{scale}}\right) - n\log(\mathrm{scale}).\]

There is no closed-form MLE; we typically optimize it numerically.

# Numerically verify the right-tail constants using SciPy

x = np.logspace(0, 6, 200)  # 1 ... 1e6
tail_pdf_const = landau.pdf(x) * x**2
tail_sf_const = landau.sf(x) * x
limit = 2.0 / np.pi

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=(r"x^2 f(x) → 2/π", r"x·P(X>x) → 2/π"),
)
fig.add_trace(go.Scatter(x=x, y=tail_pdf_const, mode="lines"), row=1, col=1)
fig.add_hline(y=limit, line=dict(dash="dash"), row=1, col=1)
fig.update_xaxes(title_text="x", type="log", row=1, col=1)
fig.update_yaxes(title_text=r"x^2 f(x)", row=1, col=1)

fig.add_trace(go.Scatter(x=x, y=tail_sf_const, mode="lines"), row=1, col=2)
fig.add_hline(y=limit, line=dict(dash="dash"), row=1, col=2)
fig.update_xaxes(title_text="x", type="log", row=1, col=2)
fig.update_yaxes(title_text=r"x·sf(x)", row=1, col=2)

fig.update_layout(width=1050, height=420, showlegend=False)
fig.show()

print("2/π ≈", limit)
print("x^2 f(x) at 1e6:", float(tail_pdf_const[-1]))
print("x·sf(x) at 1e6:", float(tail_sf_const[-1]))

2/π ≈ 0.6366197723675814
x^2 f(x) at 1e6: 0.6366302229376012
x·sf(x) at 1e6: 0.6366252002772144

def landau_loglik(x: np.ndarray, loc: float, scale: float) -> float:
    x = np.asarray(x, dtype=float)
    if scale <= 0:
        return -np.inf
    return float(np.sum(landau.logpdf(x, loc=loc, scale=scale)))


def landau_loc_scale_init(x: np.ndarray) -> tuple[float, float]:
    """Robust initializer based on median and IQR scaling."""

    x = np.asarray(x, dtype=float)
    med = float(np.median(x))
    q25, q75 = np.quantile(x, [0.25, 0.75])
    scale_init = float(max((q75 - q25) / iqr0, 1e-6))
    loc_init = float(med - scale_init * median0)
    return loc_init, scale_init


def landau_mle_scipy(x: np.ndarray) -> tuple[float, float]:
    """MLE via SciPy optimizer on (loc, log_scale)."""

    x = np.asarray(x, dtype=float)

    def nll(theta: np.ndarray) -> float:
        loc = float(theta[0])
        scale = float(np.exp(theta[1]))
        return -landau_loglik(x, loc=loc, scale=scale)

    loc0, scale0 = landau_loc_scale_init(x)

    res = optimize.minimize(
        nll,
        x0=np.array([loc0, np.log(scale0)]),
        method="Nelder-Mead",
        options={"maxiter": 4000},
    )
    loc_hat, log_scale_hat = res.x
    return float(loc_hat), float(np.exp(log_scale_hat))


# Demonstrate likelihood estimation
true_loc, true_scale = 0.8, 1.2
x_sample = landau.rvs(loc=true_loc, scale=true_scale, size=600, random_state=rng)

loc_hat, scale_hat = landau_mle_scipy(x_sample)
loc_init, scale_init = landau_loc_scale_init(x_sample)

print("true (loc, scale) =", (true_loc, true_scale))
print("init (loc, scale) =", (loc_init, scale_init))
print("MLE  (loc, scale) =", (loc_hat, scale_hat))

true (loc, scale) = (0.8, 1.2)
init (loc, scale) = (0.822319499230927, 1.2222099320256232)
MLE  (loc, scale) = (0.845424930556134, 1.1921741420882102)

7) Sampling & simulation (NumPy-only)#

The Landau distribution is a special case of a stable distribution (\(\alpha=1\)). A practical way to sample it is the Chambers–Mallows–Stuck (CMS) method.

For the standard Landau, CMS (specialized to \(\alpha=1\) and maximal right skew) can be written:

  1. Draw \(U \sim \mathrm{Unif}(-\pi/2,\pi/2)\) and \(W \sim \mathrm{Exp}(1)\) independently.

  2. Return

\[X = \frac{2}{\pi}\left[(\tfrac{\pi}{2}+U)\tan U\; -\; \log\left(\frac{(\pi/2)\,W\cos U}{\tfrac{\pi}{2}+U}\right)\right].\]

Then apply location/scale: \(\;\mathrm{loc} + \mathrm{scale}\,X\).

Numerically, we clip \(U\) away from \(\pm\pi/2\) to avoid overflow in tan and cos.

def landau_rvs_numpy(
    rng: np.random.Generator,
    size: int,
    loc: float = 0.0,
    scale: float = 1.0,
    eps: float = 1e-12,
) -> np.ndarray:
    """NumPy-only sampler for Landau via the CMS method.

    Parameters match SciPy's (loc, scale) transform.
    """

    if scale <= 0:
        raise ValueError("scale must be > 0")
    if not (0 < eps < 1e-2):
        raise ValueError("eps must be small and positive")

    u = rng.random(size)
    u = np.clip(u, eps, 1.0 - eps)
    U = (u - 0.5) * np.pi  # in (-π/2, π/2)
    W = rng.exponential(1.0, size=size)

    x0 = (2.0 / np.pi) * (
        (np.pi / 2.0 + U) * np.tan(U)
        - np.log(((np.pi / 2.0) * W * np.cos(U)) / (np.pi / 2.0 + U))
    )

    return loc + scale * x0


# Quick validation against SciPy percentiles
n = 250_000
x = landau_rvs_numpy(rng, n)
qs = np.array([0.01, 0.1, 0.5, 0.9, 0.99])
emp = np.quantile(x, qs)
the = landau.ppf(qs)

print("quantiles:", qs)
print("empirical:", emp)
print("theory   :", the)
print("diff     :", emp - the)

quantiles: [0.01 0.1  0.5  0.9  0.99]
empirical: [-1.6281 -0.9793  0.5848  7.0574 62.8577]
theory   : [-1.6275 -0.9828  0.5756  7.1287 66.0205]
diff     : [-0.0006  0.0035  0.0091 -0.0713 -3.1629]

8) Visualization#

We’ll visualize:

  • the PDF and CDF of the standard Landau

  • a histogram of Monte Carlo samples with PDF overlay

  • the instability of the sample mean (a symptom of the missing expectation)

# PDF and CDF (standard Landau)

xgrid = np.linspace(-5, 20, 5000)

fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))
fig.add_trace(go.Scatter(x=xgrid, y=landau.pdf(xgrid), mode="lines", name="pdf"), row=1, col=1)
fig.add_trace(go.Scatter(x=xgrid, y=landau.cdf(xgrid), mode="lines", name="cdf"), row=1, col=2)

fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(width=950, height=380, showlegend=False)
fig.show()

print("0.99 quantile:", float(landau.ppf(0.99)))

0.99 quantile: 66.02051286847637

# Monte Carlo: histogram + PDF overlay (clipped for readability)

n = 200_000
x = landau_rvs_numpy(rng, n)

clip_lo, clip_hi = -5.0, 20.0
x_vis = x[(x >= clip_lo) & (x <= clip_hi)]
fraction_clipped = 1.0 - (len(x_vis) / len(x))

xbins = np.linspace(clip_lo, clip_hi, 140)

fig = go.Figure()
fig.add_trace(
    go.Histogram(
        x=x_vis,
        xbins=dict(start=clip_lo, end=clip_hi, size=xbins[1] - xbins[0]),
        histnorm="probability density",
        name="samples",
    )
)
fig.add_trace(
    go.Scatter(
        x=xgrid,
        y=landau.pdf(xgrid),
        mode="lines",
        line=dict(width=3),
        name="true pdf",
    )
)
fig.update_layout(
    title=f"Landau samples (n={n:,}) — histogram clipped to [{clip_lo}, {clip_hi}]",
    xaxis_title="x",
    yaxis_title="density",
    width=950,
    height=450,
)
fig.show()

print("fraction clipped:", fraction_clipped)

fraction clipped: 0.034275000000000055

# Running mean is unstable; robust location estimates behave better

n = 60_000
x = landau_rvs_numpy(rng, n)

running_mean = np.cumsum(x) / (np.arange(n) + 1)

# For comparison: mean after clipping the extreme right tail
clip = 50.0
x_clip = np.clip(x, -np.inf, clip)
running_mean_clip = np.cumsum(x_clip) / (np.arange(n) + 1)

# Approximate running median using block medians
block = 250
m = n // block
block_medians = np.array([np.median(x[i * block : (i + 1) * block]) for i in range(m)])
avg_block_median = np.cumsum(block_medians) / (np.arange(m) + 1)

fig = make_subplots(
    rows=2,
    cols=1,
    vertical_spacing=0.12,
    subplot_titles=(
        "Running mean (unclipped vs clipped)",
        f"Average of block medians (block={block})",
    ),
)

fig.add_trace(go.Scatter(x=np.arange(n), y=running_mean, mode="lines", name="mean"), row=1, col=1)
fig.add_trace(
    go.Scatter(x=np.arange(n), y=running_mean_clip, mode="lines", name=f"mean clipped @ {clip}"),
    row=1,
    col=1,
)
fig.update_yaxes(title_text="mean", row=1, col=1)

fig.add_trace(
    go.Scatter(x=np.arange(m) * block, y=avg_block_median, mode="lines", name="avg block median"),
    row=2,
    col=1,
)
fig.update_xaxes(title_text="sample index", row=2, col=1)
fig.update_yaxes(title_text="location", row=2, col=1)

fig.update_layout(width=950, height=650)
fig.show()

print("sample median:", float(np.median(x)))
print("sample mean:", float(np.mean(x)))
print("max sample:", float(np.max(x)))

sample median: 0.5740334667297927
sample mean: 7.450544223310818
max sample: 35810.086533710426

9) SciPy integration (scipy.stats.landau)#

Key methods (no shape parameters):

  • landau.pdf(x, loc, scale)

  • landau.cdf(x, loc, scale)

  • landau.ppf(q, loc, scale)

  • landau.rvs(loc, scale, size, random_state)

  • landau.fit(data) (MLE for loc, scale)

Remember: landau.stats(moments='mv') returns nan because the mean/variance are not finite.

# Basic SciPy usage

loc, scale = 1.3, 0.7
x = np.array([-1.0, 0.0, 1.0, 5.0])

print("pdf:", landau.pdf(x, loc=loc, scale=scale))
print("cdf:", landau.cdf(x, loc=loc, scale=scale))

samples = landau.rvs(loc=loc, scale=scale, size=5, random_state=rng)
print("rvs:", samples)

# Fit (MLE) from synthetic data
n = 8_000
true_loc, true_scale = -0.5, 1.1
data = landau.rvs(loc=true_loc, scale=true_scale, size=n, random_state=rng)

loc_fit, scale_fit = landau.fit(data)  # returns (loc, scale)
loc_iqr, scale_iqr = landau_loc_scale_init(data)

print("true (loc, scale) =", (true_loc, true_scale))
print("fit  (loc, scale) =", (float(loc_fit), float(scale_fit)))
print("IQR  (loc, scale) =", (loc_iqr, scale_iqr))

pdf: [0.     0.0247 0.4054 0.0345]
cdf: [0.     0.0023 0.2469 0.866 ]
rvs: [2.306  0.3279 0.9101 8.7543 3.5786]

true (loc, scale) = (-0.5, 1.1)
fit  (loc, scale) = (-0.5066176345155946, 1.1137226858667808)
IQR  (loc, scale) = (-0.48943289656940303, 1.1127280900386476)

10) Statistical use cases#

Hypothesis testing#

Because the mean is undefined, tests based on \(\bar{X}\) (e.g. t-tests) are inappropriate.

A robust alternative is to test location using the sample median. For a continuous distribution with median \(m\) and density \(f(m)>0\):

\[\tilde{X} \approx \mathcal{N}\left(m,\ \frac{1}{4n f(m)^2}\right).\]

For \(X\sim\mathrm{Landau}(\mathrm{loc},\mathrm{scale})\):

  • median is \(m = \mathrm{loc} + \mathrm{scale}\,m_0\) where \(m_0\) is the standard median

  • density at the median is \(f(m)=f_0(m_0)/\mathrm{scale}\)

Bayesian modeling#

  • Landau likelihood can model right-skewed, heavy-tailed errors.

  • There is no conjugacy; practical inference uses MCMC or grid methods in low dimensions.

Generative modeling#

Landau is stable (index \(\alpha=1\)), so sums of independent Landau variables remain Landau, but you must be careful about parameterizations. In the \((\mu,c)\) parameterization, addition is simple:

\[X\sim(\mu_1,c_1),\ Y\sim(\mu_2,c_2)\ \Rightarrow\ X+Y\sim(\mu_1+\mu_2,\ c_1+c_2).\]

Converting back to SciPy (loc, scale) introduces the log-shift described in Section 5.

# Location test via the sample median (known scale)

def median_z_test_landau_location(
    x: np.ndarray,
    loc_null: float,
    scale: float,
) -> tuple[float, float]:
    """Approximate two-sided z-test for the location using the sample median."""

    x = np.asarray(x, dtype=float)
    if scale <= 0:
        raise ValueError("scale must be > 0")

    n = x.size
    med = float(np.median(x))

    # Under H0: median = loc_null + scale * median0
    m_null = loc_null + scale * median0
    se = scale / (2.0 * np.sqrt(n) * f0_at_median)

    z = (med - m_null) / se
    p = 2.0 * (1.0 - norm.cdf(abs(z)))
    return float(z), float(p)


n = 401
scale = 1.0

# Under H0
x_h0 = landau.rvs(loc=0.0, scale=scale, size=n, random_state=rng)
print("H0 example:", median_z_test_landau_location(x_h0, loc_null=0.0, scale=scale))

# Under H1 (shifted)
x_h1 = landau.rvs(loc=0.8, scale=scale, size=n, random_state=rng)
print("H1 example:", median_z_test_landau_location(x_h1, loc_null=0.0, scale=scale))

H0 example: (0.10679392033603453, 0.9149524707614403)
H1 example: (7.318171989669451, 2.5135449277513544e-13)

# A simple Bayesian example: posterior over loc with known scale (grid approximation)

tau = 2.0  # prior std for loc
true_loc, scale = 1.0, 1.0

x = landau.rvs(loc=true_loc, scale=scale, size=60, random_state=rng)

# Grid over loc near a robust center
center = float(np.median(x) - scale * median0)
grid = np.linspace(center - 6.0, center + 6.0, 5001)
dx = grid[1] - grid[0]

log_prior = norm.logpdf(grid, loc=0.0, scale=tau)
log_like = np.sum(landau.logpdf(x[:, None], loc=grid[None, :], scale=scale), axis=0)
log_post = log_prior + log_like

# stabilize + normalize
log_post -= np.max(log_post)
post = np.exp(log_post)
post /= np.trapz(post, grid)

post_cdf = np.cumsum(post) * dx
post_cdf /= post_cdf[-1]

loc_map = float(grid[np.argmax(post)])
loc_med = float(np.interp(0.5, post_cdf, grid))

fig = go.Figure()
fig.add_trace(go.Scatter(x=grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=true_loc, line=dict(dash="dash"), annotation_text="true loc")
fig.add_vline(x=loc_map, line=dict(dash="dot"), annotation_text="MAP")
fig.add_vline(x=loc_med, line=dict(dash="dot"), annotation_text="posterior median")
fig.update_layout(
    title="Posterior over loc (Landau likelihood, Normal prior; scale known)",
    xaxis_title="loc",
    yaxis_title="density",
    width=950,
    height=420,
)
fig.show()

print("true loc:", true_loc)
print("MAP:", loc_map)
print("posterior median:", loc_med)

true loc: 1.0
MAP: 1.2550943293575374
posterior median: 1.2440297824189852

# Generative property: (μ, c) stable addition vs SciPy's (loc, scale)

def loc_to_mu(loc: float, scale: float) -> float:
    if scale <= 0:
        raise ValueError("scale must be > 0")
    return float(loc - (2.0 * scale / np.pi) * np.log(scale))


def mu_to_loc(mu: float, scale: float) -> float:
    if scale <= 0:
        raise ValueError("scale must be > 0")
    return float(mu + (2.0 * scale / np.pi) * np.log(scale))


loc1, c1 = 0.0, 1.0
loc2, c2 = 0.0, 2.0

n = 250_000
x = landau_rvs_numpy(rng, n, loc=loc1, scale=c1)
y = landau_rvs_numpy(rng, n, loc=loc2, scale=c2)
z = x + y

# Convert to (μ, c), add, then convert back
mu1, mu2 = loc_to_mu(loc1, c1), loc_to_mu(loc2, c2)
mu_z = mu1 + mu2
c_z = c1 + c2
loc_z = mu_to_loc(mu_z, c_z)

qs = [0.1, 0.5, 0.9]
q_emp = np.quantile(z, qs)
q_naive = landau.ppf(qs, loc=loc1 + loc2, scale=c_z)
q_corr = landau.ppf(qs, loc=loc_z, scale=c_z)

print("empirical quantiles:", q_emp)
print("naive (loc add)  :", q_naive)
print("corrected        :", q_corr)
print("corrected loc_z  :", loc_z)

empirical quantiles: [-1.74    2.9265 22.2747]
naive (loc add)  : [-2.9485  1.7269 21.386 ]
corrected        : [-1.7329  2.9425 22.6017]
corrected loc_z  : 1.2156525147857526

11) Pitfalls#

  • Do not use the sample mean/variance as estimators; they are unstable because the corresponding moments do not exist.

  • Underflow in the far left tail: pdf(x) can become numerically 0 for large negative x. Prefer logpdf for inference.

  • Parameter validity: scale must be strictly positive.

  • Fitting: MLE can be sensitive with small samples because a few large values dominate the likelihood. Use robust initializations (median/IQR).

  • Visualization: histograms need clipping or log axes to avoid being dominated by rare extreme values.

12) Summary#

  • landau is a continuous, highly right-skewed distribution with a power-law right tail.

  • Mean and variance do not exist; prefer median/quantiles and likelihood-based modeling.

  • Sampling is convenient via the CMS stable-variable algorithm (NumPy-only).

  • SciPy’s scipy.stats.landau provides accurate pdf/cdf/ppf/rvs/fit implementations.

  • When adding Landau variables, be careful about parameterization: the \((\mu,c)\) form adds cleanly, and it maps to SciPy’s (loc, scale) with a log shift.

Contents